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Abstract 

Schemes with the second-order approximation in time are considered for 
numerical solving the Cauchy problem for an evolutionary equation of first or¬ 
der with a self-adjoint operator. The implicit two-level scheme based on the 
Fade polynomial approximation is unconditionally stable. It demonstrates good 
asymptotic properties in time and provides an adequate evolution in time for 
individual harmonics of the solution (has spectral mimetic stability). In fact, 
the only drawback of this scheme is the necessity to solve an equation with an 
operator polynomial of second degree at each time level. We consider modifi¬ 
cations of these schemes, which are based on solving equations with operator 
polynomials of first degree. Such computational implementations occur, for 
example, if we apply the fully implicit two-level scheme (the backward Euler 
scheme). A three-level modification of the SM-stable scheme is proposed. Its 
unconditional stability is established in the corresponding norms. The emphasis 
is on the scheme, where the numerical algorithm involves two stages, namely, the 
backward Euler scheme of first order at the first (prediction) stage and the fol¬ 
lowing correction of the approximate solution using a factorized operator. The 
SM-stability is established for the proposed scheme. To illustrate the theoretical 
results of the work, a model problem is solved numerically. 

Keywords: Evolutionary equation of first order, the Cauchy problem, finite 
difference schemes. Fade approximation, SM-stability, factorized scheme 
PACS: 02.30.Jr, 02.60.Lj, 02.70.Bf 
2000 MSC: 65J08, 65M06, 65M12 


’^This work was supported by the Russian Foundation for Basic Research (project 14-01- 
00785) 

* Correspondibg author. 

Email address: vabishchevichSgmail. com (P.N. Vabishchevich) 


Preprint submitted to Elsevier 


April 17, 2015 



1. Introduction 


In the practice of scientific computations for numerical solving transient 
boundary value problems, focus is on computational algorithms of higher accu¬ 
racy (see, e.g., ^[5]). Along with increasing the accuracy of approximation in 
space, researchers try to improve the accuracy of approximation in time using 
primarily numerical methods developed for ordinary differential equations El El- 
Concerning unsteady problems governed by partial differential equations, we are 
interested in numerical methods to solve the Cauchy problem for stiff systems 
of ordinary differential equations laiiiii. 

Increasing of the accuracy of the approximate solution to time-dependent 
problems is achieved in different ways. For two-level schemes, where the solu¬ 
tion at two consecutive time levels is involved, polynomial approximations are 
explicitly or implicitly employed for operators of difference schemes. The Runge- 
Kutta methods 011 ] are well-known examples of such schemes widely used in 
modern computing practice. The main feature of multilevel schemes (multistep 
methods) results in the approximation of time derivatives with a higher accu¬ 
racy on a multi-point stencil. Multistep methods based on numerical backward 
differentiation formulas [5] should be mentioned as a typical example. 

In the transition from continuous problems to discrete ones, we seek to in¬ 
herit the most important properties of the differential boundary value problem 
under the consideration. For time-dependent problems, emphasis is on the main 
criterion to well-posedness of the grid problem, i.e., the stability of the differ¬ 
ence solution with respect to small perturbations in the initial data, boundary 
conditions, right-hand sides and coefficients of equations uni. 

Various classes of stable difference schemes can be used. Among stable differ¬ 
ence schemes, we search for a scheme that is optimal in sense of some additional 
criteria. In the theory of difference schemes, the class of asymptotically stable 
difference schemes [Hill] was highlighted as the class of methods ensuring the 
correct behavior of the approximate solution at long time. In the theory of 
numerical methods for solving systems of ordinary differential equations 00, 
the concept of L-stable methods was introduced, where an asymptotic behavior 
of the approximate solution is also investigated at long time. 

Spectral Mimetic (SM) properties of operator-difference schemes for numer¬ 
ical solving the Cauchy problem for evolutionary equations of first order are 
associated with the time-evolution of individual harmonics of the solution. A 
study of spectral characteristics allows to select more acceptable approximations 
in time. The concept of SM-stable difference schemes has been introduced in 
m for two-level schemes of higher accuracy for solving the Cauchy problem for 
evolutionary first-order equations. The long-time behavior of individual har¬ 
monic of the approximate solution was studied for problems with a self-adjoint 
operator using various implicit schemes, which are based on the Fade approx¬ 
imations |14j with the corresponding operator exponential. In [laiii] , it was 
shown that the best way to obtain SM-properties for such problems is to apply 
polynomial approximations. 

Among two-level implicit schemes of higher accuracy order for evolutionary 
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equations with skew-symmetric operators that were obtained on the basis of the 
Fade approximations, symmetric schemes m demonstrate some advantages. 
For evolutionary equations with general operators, it is possible to construct 
additive schemes (splitting schemes) via splitting the problem operator into 
self-adjoint and skew-symmetric components m- 

The analysis of schemes for problems with self-adjoint operators conducted 
in the mentioned above papers allowed to highlight the class of SM-stable dif¬ 
ference schemes based on the Fade polynomial approximation. The implicit 
two-level scheme of second order approximation in time constructed using the 
Fade polyrromial approximation is unconditionally SM-stable. The computa¬ 
tional implementation of this scheme involves the solution of an equation with 
an operator polynomial of second degree at each time step. It seems reason¬ 
able to modify these schemes in order to employ more simple equations with an 
operator polynomial of first degree. For example, we can use the fully implicit 
(backward Euler) two-level scheme. In this case, at each time level, we solve the 
problem with a factorized operator, where each factor is an operator polynomial 
of first degree. When considering systems of ordinary differential equations, in 
this case, we speak of diagonally implicit Runge-Kutta (DIRK methods) [311]. 
Factorized SM-stable schemes were investigated in the work [18]. It was found 
that among the schemes with maximal accuracy, in fact, there is no SM-stable 
scheme. 

In the present paper, we construct factorized schemes of second-order ap¬ 
proximation in time to the Cauchy problem for the evolutionary equation of 
first order with a self-adjoint operator. They are based on a modification of the 
implicit two-level scheme based on the Fade polynomial approximation, which 
has the second-order approximation in time, is unconditionally stable, demon¬ 
strates good asymptotic properties in time and adequate long-time behavior of 
individual harmonics of the solution (spectral mimetic stability). 

In the first modification, instead of two-level scheme, we employ three-level 
scheme. Here a part of the approximation for the time derivative is taken 
from the previous time level. A similar idea we already used (see, e.g., [T5j) to 
increase the accuracy of explicit-implicit schemes. The unconditional stability of 
the three-level scheme is shown in the corresponding norms. A more interesting 
modification employs the two-stage calculations, namely, the backward Euler 
scheme of first order at the prediction stage and the following correction of the 
approximate solution via solving the problem with a factorized operator. This 
scheme is unconditionally SM-stable. 

The paper is organized as follows. In Section 2, we consider unconditionally 
SM-stable scheme of second-order approximation in time for solving the Cauchy 
problem for the evolutionary equation of first order. The factorized scheme 
based on the three-level modification is investigated in Section 3. Section 4 is 
the core of our work. Here we study the SM-stable factorized scheme. Numerical 
experiments for a model problem with a one-dimensional parabolic equation are 
presented in Section 5. 
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2. SM-stable scheme of second order 


Let us consider a finite-dimensional real Hilbert space H, where the scalar 
product and the norm are (•, •) and || • ||, respectively. Suppose that u{t) (0 < t < 
T,T > 0) is defined as the solution of the Cauchy problem for the homogeneous 
evolutionary first-order equation: 

du 

——\-Au = 0, 0 < t < T, (1) 

at 

u(0) = u°. (2) 

In equation Q, H is linear and, for simplicity, independent of t (constant) 
operator acting from H into H {Ad/dt = d/dtA). In addition, assume that the 
operator is self-adjoint positive definite: 


A = A*>dE, S>0, 


( 3 ) 


where E is the identity operator. 

For problem 0-@> the following elementary estimate for stability of the 
solution with respect to the initial data holds: 

||M(t)|| < exp(-5t)||w°||. (4) 

SM-properties for approximations in time are connected [18] with inheriting the 
basic properties of the differential problem. 

To solve numerically the problem 0, 0, we employ two-level difference 
schemes. Define the uniform grid in time with a step r as 

Wt = Wt U {T} = {f" = nr, n = 0,1,..., N, tN = T} 

and denote y" = y{tA), where = nr. We want to pass from a time level t" to 
the next time level The exact solution has the form 

n(a:, = exp(—rH)'u(a::, t"). (5) 

The two-level scheme corresponding to 0 may be written in the form 

= n = 0,l,..., (6) 

where S is the transition operator from the current time level to the next one. 
In the general case, the operator S may depend on n. We restrict the considera¬ 
tion to the difference approximations with respect to time for the homogeneous 
problem 0,0 that lead to the transition operator 

S = s(rH), (7) 

where s(z) is called a stability function. 
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The use of additional information about properties of a problem operator 
allows to specify requirements for the choice of approximations in time. Taking 
into account the inequality ([^, we get a more accurate stability estimate 

||u(a;,r+^)|| < exp(-(5r)||u(x,t")||. 

The inheritance of this property of the problem means that the following 

estimate must be true: 


< p||2/”||, ||s(rA)||<p, p = exp(-(5T). (8) 


In this case, we say m that this operator-difference (semi-discrete) scheme is 
p-stable. 

The properties discussed above concern the bounds of the spectrum of the 
transition operator. We want also to study some other qualitative indicators of 
the approximate solution behavior. In solving unsteady problems, we focuses 
on the long-time asymptotic behavior of the solution. For the model problem 
under examination, the solution decreases to zero as t —>■ oo. If the approximate 
solution preserves this property, we say that the approximate solution is asymp¬ 
totically stable (with respect to time). In the case of the Cauchy problem for 
systems of ODEs, this property of the approximate solution is called L-stability. 
If 


lim s{9) = 0, (9) 

6—>oo 


then the stable (p-stable) difference scheme Q is said to be asymptotically 
stable. 

For the linear problem 0.0 with a self-adjoint positive definite operator 
A, the solution may be written as the superposition of individual harmonics, 
which are associated with their own eigenvalues. In m, the choice of approx¬ 
imations in time is subjected to the requirement of the appropriate behavior 
in time of all points of the spectrum. In this case, we introduce SM-property 
(Spectral Mimetic) of schemes, i.e., such a scheme is said to be SM-stable. A 
difference scheme is called SM-stable if it is p-stable and asymptotically stable. 
The additional requirement is the spectral monotonicity, namely, the function 
s{0) is monotonically decreasing. This means that the harmonics with higher 
indexes decay more rapidly than the harmonics with lower indexes. 

Two-level schemes of higher-order approximations for time-dependent linear 
problems can be conveniently constructed on the basis of the Fade approxi¬ 
mations of the operator (matrix) exponential function exp(—rA). Other ap¬ 
proaches are discussed in EOllsi]. In the case of nonlinear systems of ODEs, 
such approximations correspond to various variants of Runge-Kutta methods 

[DiaiT]. 

The Fade approximation of the function exp(—z) is 


exp(-z) = Rim{z) + £>( 2 '+”"+^), 




Plmjz) 
Qlm {z) 


( 10 ) 


5 



where Pim{z) and Qim(z) are polynomials of degrees I and m, respectively. 
These polynomials have mi the form 




ll 

{I + m)\ 


{l + m-ky 

h ’ 


Qlm (^) 


ml 

{I + m)\ 


(^ + ™ ~ ^)- ,fc 
kl{m — k)l 


Only the schemes with I = 0, i.e., 


Pomi,z) — / \ ’ Pom{z) — 1 

^Om \^) 

are SM-stable. In this case, the two-level difference scheme for the problem Q, 
([^ is 

n+l _ n I 

QomirA)- -h-(Qom(TA) - S)?/" = 0, n = 0, 1,..., (11) 

T T 

where the function 

m ^ 

Q«r,M = E 

k^O 

is a truncated Taylor series for exp(z). Figure presents the corresponding 
stability functions for m = 1,2,3. 

When considering the Cauchy problem with an inhomogeneous right-hand 
side, instead of Q, we solve the equation 

^+Au = fit), 0<t<T, (12) 

The SM-stable scheme of hrst order (m = 1 in 0) may be written in the form 

n+l _ n 

{E -t tA) ^^ -h Ay" = yj". (13) 

T 

Taking into account that the exact solution of (i), (HI) satisfies the representa¬ 
tion 


i(t"+i) = exp 1^—T ^ Aaj u{E) 

+ It exp f-(t"+^-6>)^ A„'j/(6l)d6», 




for the right-hand side of equations (13), we can take = /(t"^^^^). 
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Figure 1: SM-stable Fade approximations 


The object of our study is the factorized scheme of second order (m = 2 in 
0) , where 


1 \ 7 /"'+! — 7 /*^ 

E + TA+-T^A^^y y 


+ [A+-TA^]y- = ^-. 


(14) 


For the right-hand side, it seems reasonable to use the approximation 


= + /(r+i/2). 


The main drawback of the scheme (14) is associated with solving the problem 

1 


E + tA+ 2/”+^ = X” 

at each time level, where we have A^. The computational implementation of 
the backward Euler scheme (131 involves the solution of the essentially simpler 
problem (E + = x"- A fundamental simplification in solving the 

problem at a new time level is connected with using the factorized scheme, 
where we solve the problem 


(E + aiTA)(E + cr2T74)y"+^ = y” 
with positive constants aa, a = 1, 2. 


(15) 
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3. Three-level factorized scheme 


To construct factorized schemes, we restrict our study to the case with the 
identical weighting parameters in (151, i.e., the solution at the new time level is 
evaluated from the equation 


{E + aTAfy^+^ =X^. 

Let us rewrite the basic scheme ( |l4| in the form 

n+l _ / 1 

{E + arAf' - \- i E + tA H— A? — (E -\- arA) 

T \ 2 

+ {A+^rAAy- = y,\ 


2 \ 


(16) 


Factorized schemes are based on some modification of the second term in the 
left-hand side of (16). In view of 


E + tA+ -{E + utAY = (1 - 2a)TA + ~ ^ ^ ’ 

to preserve the second-order approximation, it is sufficient to approximate the 
difference derivative — y")/!" with the first order with respect to r. 

Consider the case of modification, where 


y"+i - y" 


y - y 


n—1 


Such a transition from a two-level scheme to a three-level scheme was used in m 
in the construction of explicit schemes for parabolic and hyperbolic equations. 


Instead of (14), we apply the scheme 

,,"-1-1 _/ 1 

{E + arAY - \- [ E + tA -\—— {E + arA) 

T \ 2 

+ U+ItaA y- = ip\ 


Taking into account that 

yn+l_yn yn+l_yn-l ^ yU+l _ 2yn yU-l 


2 \ 


(17) 


yn_yn-l yn+l_yn-l ^ yU+l _ 2yn yU-l 

T 2 t 2 ’’ 

rewrite it in the canonical form m for three-level operator-difference schemes: 


n+l _ n ,,"-1-1 _ 2?/” -I- ~ 

- y_ ^ j^y -^ 


( 18 ) 




















For the operators A, B,D, we have 


A = A+-tA\ 


B = E + tA + -t^A^, 


D = -{E + {‘ia- 1 )tA + 20-2 _ 


1 




The necessary and sufficient conditions for the stability of the three-level 
scheme (18) with constant self-adjoint operators A,B,D, where B > 0 and 
A > 0 , are formulated (see (TUI [TTl [52] ) in the form of the inequality 

^2 


D > —A. 
4 


(19) 


We have 


G = D- —A = - {E+[Aa- - \TA+[2a^ - - \ t^A^ . 


The inequality (19) holds for a > y^3/8. 

Theorem 1. The three-level factorized scheme Jl?| ) is unconditionally stable 
under the restriction a > ^^3/8. The following a priori estimate holds: 


£n+l < ^n + 2 


( 20 ) 


where 




y"+i -h y" 

2 

y"-Hl _ y" 

- 1 - 

A 

2 

r 


Proof. By 


y" = + 2 y” + 2 /”-i) - - 2 y" + 


we can rewrite (18) as 

n-l-l _ „ra-l _ 97," 4- 7,"-l .^7,"-l-l _ 97," 4_ 7,"-l 

- y —-fy+y—^^y- ^jL_±y_—^ n_ (21) 

2 r 4 ^ ^ 


Let 


1 _ ,,"-1 

^" = ;5(y" + y”-'), = ^ 

Z T 


then ( 21 ) can be written in the form 

yjn+l_^n I 


B 


yjU+l _|_ yjU 


+ G- 


+ iA(?;”+i +z;”) =(^". 


( 22 ) 
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Multiplying scalarly both sides of \22\ by 

2 (u"+i - u”) = + u;"), 

we get the equality 

+ w’"), u;”+i + w”) + (G(u>”+i - w"), + w’") 

+ (I(u"+i + u"), u”+i - u”) = t((^", u;"+i + w^). 

To estimate the right-hand side, we use the inequality 

u;”+i + u;") < ^(S(u;"+i + u;”),u;”+i + u.") + V",O- 


( 23 ) 


This makes it possible to get from (23) the inequality 


where we use the notation 

£„ = (Iu’^,u") + (Gu;",u;”). 


(24) 


The inequality (24) is the desired a priori estimate (20). 


4. SM-stable factorized schemes 


The second possibility to modify the scheme (16) is connected with the 
solution of an auxiliary problem. In this case 


_ yU 


in calculating the auxiliary quantity which approximates . This stage 
of predicting the solution at the new time level can be naturally implemented 


on the basis of SM-stable scheme of first order. According to (131, we put 


~n+l _ n 

{E + tAY- -^ + Ay^ = 


(25) 


where, for example, = /(t"“''^/^). At the correction stage, we apply the 
factorized scheme 


n-l-l _ ,,n / 1 \ 

{E + arAf- - — + {e + tA+^YA^ - {E + aTAf\ - - — 

+ fA+^-TAAy- = ^\ 


(26) 


Let us formulate the conditions for SM-stability of the factorized scheme (25) 


( 26 ). 
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Theorem 2. The factorized scheme (25), (26) is unconditionally SM-stable 
under the restriction a > If \/2• 


Proof. From (25), we have 


y 


y 


= -{{E + TA)-\Ay^-tp^). 


Substitution in (26) results in 

n+1 _ , 


(E + tA){E + arA) 


2y 


+ {E + tA) (^E + Ay” = 


— ( (1 — 2cr)rA + 

) 


2 


r^A^ Ay 


(27) 


The two-level scheme (27) may be written as 

y”+l = ^y” + rV’". 

For the stability function s(z), we have the representation 

1 -I- 2az + (ct^ — 0.5)z^ 


s{z) = 


1 -I- (1 -I- 2a)z + (cr^ -I- 2a)z'^ + a'^z^ 


(28) 


It is easy to see that the condition of asymptotic stability (|^ holds for (28). The 
most important condition of SM-stability that is associated with monotonous 
decreasing the function s{z) with z > 0 is satisfied for cr > l/v^. This proves 
the theorem. 

The limiting value a = l/v^ provide the best approximation properties. It 
can be observed in Figure]^ Note that the SM-stable scheme (251, (26) is based 
on the factorization into three operators of a simpler structure, whereas in the 
three-level scheme 0. we have only two factors. The main advantage of the 
scheme (25), (p^ is in its absolute SM-stability. 


5. Numerical experiments 

Now we will illustrate the possibilities of the above factorized schemes of 
second-order accuracy in time on predictions of a model parabolic problem. 
Consider a boundary value problem for the equation 


du d'^u 

m~dx^ 


= f{x, t), 0 < X < I, 0 < t < T. 


(29) 


The equation (29) is supplemented with the boundary and initial conditions: 


u(0, t) = 0, u(I, t) = 0, 0 < t < T, 

u{x, 0) = 0, 0 < X < I. 


(30) 

(31) 
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Figure 2: SM-stable factorized scheme 


The problem (29)-(311 is considered with T = 0.5, when the exact solution is 
u{x, t) = a;(l — x)(l — exp(—5t)). 


To solve numerically the problem (29)-(31), we introduce a difference ap¬ 
proximation in space. Within the unit interval, we introduce a uniform grid 
with step h: 

uj = coLi du} = {x \ X = Xi = ih, i = 0,1,..., M, Mh = 1}, 

where oj is the set of internal nodes {i = — 1) of the grid. For 

u{x,t), X G uj, we have the equation 

^ + Ay = f{x,t), xGuj. 

Here the discrete operator A is defined as follows: 

' 2v{x) — v{x + h) 


/l2 

^ 2 u(a:) - v{x - h) - v{x + h) 

2v{x) — v{x — h) 
h? ’ 


X = Xi, 


X = Xi, i = 2, 3,..., M — 2, 


X = XM-1- 
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The accuracy of the approximate solution is estimated in the grid norm of the 
space L 2 {ijj), where the scalar product and the norm are defined as 

{v,w) = ''^^v{x)w(x)h, ||u|| = 

xGoj 


Thus, for the error of the approximate solution y" — u(x,t^), x G ui, we study 

eW = b” - t = 

In our numerical experiments, emphasis is on the accuracy in time of the 
difference schemes under consideration. Because of this, we present predictions 
conducted using the same spatial grid M = 10. The calculations are performed 
using different grids in time, namely, N = 10,20,40. 

The accuracy of conventional two-level difference schemes is presented in 
Figure and Figure Figure demonstrates the results for the backward 
Euler scheme (13) (the scheme of first-order accuracy). Similar numerical data 
for the Crank-Nicolson scheme 


1 A - 


+ Ay^ = 


are depicted in Figure]^ (the scheme of second-order accuracy). 



The main object of our study is the modihed SM-stable scheme of second- 
order accuracy (14). The accuracy of this scheme is shown in Figure]^ The ap¬ 
proximate solution of the model problem (29)-(311 obtained via the scheme (141 
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Figure 4: Crank-Nicolson scheme 
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Figure 5: SM-stable scheme of second-order accuracy 































demonstrates higher accuracy in comparison with the Crank-Nicolson scheme 
(compare Figurej^and Figure]^. 

When using the three-level factorized scheme 0, we must preliminarily cal¬ 
culate the approximate solution at the first time level y^. Figurej^demonstrates 
accuracy of predictions, where the exact is used. In practical calculations, 
we should focus on the schemes of the second-order accuracy for evaluating y^. 
Thus, it seems natural to use the Crank-Nicolson scheme. In this case, the 
accuracy of calculations for the model problem (29|-(31) is shown in Figure]^ 



Figure 6: Three-level factorized scheme; exact 


Similar results obtained using the SM-stable factorized scheme (25), (261 are 
presented in Figure]^ and Figure]^ The optimum value (cr = l/v% Figure]^ 
provides significantly higher accuracy in comparison with the case, where we 
employ the scheme with the twice value of the parameter a (Figure]^. 

These results demonstrate the efficiency of factorized variants of the SM- 
stable scheme with the second-order accuracy for solving parabolic boundary 
value problems. 
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Figure 9: SM-stable factorized scheme: cr = \/2 
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